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SOME  ERROR  PROBABILITIES  FOR  THE  ASSOCIATION  OF 
PASSIVE  DF  MEASUREMENTS  WITH  RADAR  RETURNS 

INTRODUCTION 

A  great  deal  of  attention  is  currently  focused  on  multiple-sensor  surveillance  systems. 
One  important  aspect  of  such  systems  is  the  correct  association  of  data  from  different 
sensors  with  the  same  target.  Of  particular  interest  here  is  the  problem  of  associating  a  set  of 
measurements  from  a  passive  direction-finding  (DF)  system  with  the  correct  target  in  a 
surveillance  radar  system.  This  can  be  viewed  as  a  problem  in  statistical  pattern  recognition 
in  which  a  set  of  DF  measurements  is  to  be  assigned  to  a  class  whose  characteristics  are 
determined  from  measurements  of  a  radar  target.  This  classification  is  performed  by 
choosing  the  radar  target  that  will  minimize  a  discriminant  which  is  a  function  of  the  DF 
data  and  the  data  on  the  radar  target. 

One  difficulty  inherent  in  this  problem  is  due  to  the  asynchronous  operation  of  the 
radar  and  DF  systems.  Because  the  targets  may  be  in  motion,  it  is  necessary  to  smooth  or 
extrapolate  in  time  the  DF  measurements,  the  radar  data,  or  both  in  order  that  informa¬ 
tion  from  both  sensors  refer  to  the  same  set  of  target  positions.  It  is  not  immediately 
obvious,  however,  which  set  of  data  should  be  smoothed.  Gerlach  [1]  assumed  simultaneous 
measurements  and  used  no  smoothing.  Coleman  [2]  smoothed  the  radar  data  only.  Bath 
[3]  smoothed  both  the  radar  and  DF  data  to  some  extent.  No  claim  was  made  by  any  of 
these  authors  that  his  choice  was  optimum,  and  in  fact  the  question  of  which  data  to 
smooth  was  not  explicitly  addressed.  The  purpose  of  this  report  is  to  bound  the  perfor¬ 
mance  obtainable  under  each  of  the  four  possible  smoothing  options  listed  in  Table  1.  In  the 
interest  of  conciseness  the  mnemonic  abbreviations  shown  in  Table  1  for  these  options  will 
be  used  freely.  Using  some  simplifying  and  unifying  assumptions  outlined  subsequently, 
exact  expressions  for  the  probability  of  classification  error  are  obtained  for  three  of  the 
options:  RDS,  RS,  and  NS.  The  probability  of  error  with  DS  is  evaluated  with  an 
importance-sampling  simulation.  Due  to  the  nature  of  the  assumptions,  these  error  probabil- 


Table  1  —  Mnemonic  Abbreviations  for 


the  Smoothing  Options 


Mnemonic 

Smoothing  Option 

RDS 

Radar  and  DF  smoothing 

RS 

Radar  smoothing  only 

DS 

DF  smoothing  only 

NS 

No  smoothing 
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ities  represent  lower  bounds  on  realistically  obtainable  performance.  These  results  for  the 
RDS  and  RS  options  have  been  described  previously  [2]  but  are  reviewed  here  in  somewhat 
more  detail.  The  results  for  the  remaining  options  are  new. 


ASSUMPTIONS 

A  number  of  assumptions  about  the  problem  are  made.  Those  based  on  earlier  work 
[2]  include  the  following: 

•  Starting  with  a  set  of  passive  measurements  and  asking  which  set  of  radar  detections 
it  should  be  paired  with.  This  assumption  assures  that  each  DF  target  is  associated  with  only 
one  radar  target  (with  the  implicit  requirement  that  radar  data  on  the  appropriate  target  be 
available).  Several  DF  targets  are  allowed  to  be  paired  with  a  single  radar  target;  for  ex¬ 
ample,  paired  with  an  aircraft  target  may  be  separate  sets  of  DF  measurements  to  each  of 
several  transmitters  on  board. 

•  Associating  an  entire  passive  measurement  set  as  a  unit  (rather  than  measurement  by 
measurement).  This  assumption  implies  an  ability  to  correctly  associate  a  new  DF  measure¬ 
ment  with  existing  DF  measurements  of  the  same  target.  It  should  allow  more  accurate  DF- 
to-radar  association  than  if  the  association  proceeded  on  a  single  DF  measurement. 

•  Basing  the  association  decision  on  minimizing  a  sum  of  squared  passive-to -radar 
bearing  differences.  This  assumption  is  based  on  a  generalized  likelihood  classifier  for 
Gaussian  random  variables  [4]  and  has  been  implicit  in  all  of  the  radar/DF  association 
schemes  known  to  this  author. 


Further  assumptions  are  made  to  simplify  the  problem.  There  are  two  targets,  and 
both  are  stationary.  This  allows  classification  without  smoothing  to  be  a  viable  option. 
Target  1  is  under  both  radar  and  DF  observation,  while  target  2  is  observed  only  by  radar. 
Due  to  the  symmetry  in  the  problem  this  assumption  can  be  made  without  loss  of  gener¬ 
ality.  Radar  and  DF  equipments  are  colocated.  Therefore  the  only  radar  measurement  of 
interest  is  azimuth.  The  targets  are  separated  in  azimuth  by  na  (defined  as  the  azimuth  of 
target  2  minus  the  azimuth  of  target  1).  The  set  of  DF  measurements  and  the  two  sets  of 
radar  measurements  contain  n  measurements  each.  All  measurements  are  unbiased  and 
independent.  The  measurement  errors  are  Gaussian,  with  the  DF  errors  having  variance 
OjjF  and  the  radar  errors  having  variance  .  Since  there  are  only  two  targets,  the  decision 
procedure  is  to  evaluate  the  discriminant  function  for  each,  form  the  difference,  and 
compare  to  zero.  Specifically,  the  procedure  is  to  compute 


(1) 


and  associate  the  DF  measurements  with  target  1  if  d  >  0  (correct  decision)  and  with 
target  2  if  d  <  0  (incorrect  decision).  In  this  equation  0DFi  is  the  ith  DF-based  azimuth 
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estimate.  If  DF  smoothing  is  not  used  (RS  and  NS),  0DFi  is  just  0DFi,  the  tth  DF  measure¬ 
ment.  If  DF  smoothing  is  used  (RDS  and  DS), 


ffDF  i  =  0DF  =  n 


n 

y,  &DFk< 

fe=i 


the  average  of  the  DF  measurements.  Similarly  is  the  ith  radar-based  estimate  of  the 
azimuth  of  target  1.  With  no  radar  smoothing  6  j,-  is  flj,-,  the  ith  radar  azimuth  measure¬ 
ment  of  target  1.  With  radar  smoothing 


fe=i 

the  average  of  the  radar  azimuth  measurements  of  target  1.  Similar  relationships  hold 
between  02l-,  6 2,  02l-  and  the  smoothing  options. 

Because  the  sign  of  d  in  Eq.  (1)  is  unaffected  by  a  scale  change  in  azimuth,  the  means 
and  standard  deviations  of  the  densities  are  normalized  by  oDF  throughout  the  subsequent 
development  except  where  noted.  The  notation 

P~Pa/o  DF 


and 


0  =  °R^a  DF 


will  be  used. 

In  the  next  four  sections  the  probability  of  classification  error  Pe  is  evaluated  for  the 
four  smoothing  options  under  the  assumptions  just  outlined. 


Pe  WITH  RADAR  AND  DF  SMOOTHING  (RDS) 

When  both  the  radar  data  and  the  DF  data  are  smoothed,  the  estimates  0DF|,  6 lf,  and 
02l-  are  all  independent  of  i.  Equation  (1)  can  therefore  be  rewritten  as 

d  =  (0DF-^2)2-(®DF-^)2 

2  -  2  (?20DF  +  <?22  “  0j2  (2) 

=  <*2-*i>K®a  +  »i>-  2^dfI- 
3 


i 
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Prom  independence  of  the  0DFi,  0li,  and  d2i  it  follows  that  the  ( 62i  -  dli)  and  (d2i  +  0li) 
are  independent.  Consequently,  the  two  factors  in  Eq.  (2)  are  independent,  and  the  product 
d  is  less  than  zero  if  and  only  if  its  two  independent  Gaussian  factors  are  of  opposite  sign. 
The  probability  of  error  can  therefore  be  expressed  as 


Pe  =  Prob  (d  <  0)  =  Prob  (6  2  -  <  0)  Prob  [(02  +  0j)  -  26 DF  >  o] 

+  Prob  (§2  -  ex  >0)Prob  [(0  2  +  6  j )  -  20Dp<o]. 

Or,  if  4>(x)  is  defined  as  the  cumulative  distribution  function  for  a  zero-mean  unit-variance 
Gaussian  random  variable. 


+ 


This  result  can  be  more  concisely  expressed  by  using  (from  Ref.  5) 


where 


erfc  (x)  =  l-  erf  (x) , 


to  obtain 
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This  function  is  plotted  in  Fig.  1  versus  >iy/n  with  a  as  a  parameter. 

Several  observations  can  be  made  regarding  this  result.  First,  the  curve  for  a  =  0  reflects 
performance  with  perfect  radar  measurements  and,  in  effect,  is  obtained  using  true  positions 
instead  of  estimates  for  the  9 1  •  and  ^2 ;  in  Eq.  (1).  Therefore  this  curve  gives  the  perfor¬ 
mance  of  an  idealized  maximum-liklihood  detector.  As  o  increases,  the  probability  of  error 
increases  slowly  until  a  approaches  unity.  This  suggests  that  the  accuracy  of  the  radar  (o^ ) 
is  not  important  as  long  as  it  is  significantly  better  than  the  accuracy  of  the  DF  equipment 

Further  insight  can  be  gained  into  the  behavior  of  this  function  by  using  a  different 
normalization.  After  substituting  the  definitions 


U  =  - 

aDF 

and 


a 


°DF 


Fig.  1  —  Performance  with  radar  and  DF  smoothing  (RDS) 
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into  Eq.  (3),  the  probability  of  error  becomes 


after  some  manipulation*.  Figure  2  shows  contours  of  constant  log10  Pe  plotted  versus 
logio  loDF/(Ma\/n)l  and  log10  \aR/((Jtay/n)\.^  It  can  be  thought  of  as  a  log-log-log  plot  of  Pe 
versus  the  two  standard  deviations  and  aR  witn  |j ua  | s/n  fixed  to  unity.  It  is  apparent 
from  this  plot  that  not  only  is  aR  unimportant  when  oR  «  oDF  but  oDF  is  unimportant 
when  o  R . 


The  question  of  which  variance  is  more  important  can  be  considered  more  rigorously 
using  sensitivities.  The  sensitivity  of  the  function  Pe  with  respect  to  aR  is  defined  as 


P  °R  dPe 
S°R  Pe  d0R 

Similarly 

P  aDF  dPe 

S  e  = - • 

P  3oni? 
DF  e  Dr 


♦As  (from  Ref.  5)  erf  ( ~x )  =  ~erf  ( x ),  Eq.  (3)  shows  Pf  to  be  an  even  function  of  p.  In  the  manipulation  to 
arrive  at  Eq.  (4)  the  sign  of  fi  is  lost  as  it  is  brought  under  the  radical  and  the  symmetry  of  Pe  with  respect 
to  fl  is  destroyed.  This  situation  was  corrected  by  replacing  po  by  |pQ|  in  Eq.  (4). 
tFigure  2  was  produced  with  the  aid  of  MACSYMA,  a  large  symbolic  manipulation  program  developed  at 
the  MIT  Laboratory  for  Computer  Science  and  supported  by  the  National  Aeronautics  and  Space  Admin¬ 
istration  under  grant  MSG  1323,  by  the  Office  of  Naval  Research  under  grant  N00014-77-C-0841,  by  the 
U.S.  Department  of  Energy  under  grant  ET-78-C-02-4687,  and  by  the  U.S.  Air  Force  under  grant 
F49620-79-C-020. 
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It  is  not  difficult  to  show  that  these  can  be  defined  equivalently  as 

Pe  d(lQgiope) 

°R  3(log10oft) 

p  W°e10Pe) 

§  e  —  - - 

°DF  9(l0gjg  Ujjp. 

These  two  quantities  are  equal  in  Fig.  2  wherever  the  contours  are  at  45°  with  respect  to 
the  axis.  The  locus  of  such  points  is  shown  in  Fig.  2  and  represents  the  boundary  between 
the  region  where  a  small  percentage  change  in  oR  is  more  significant  than  the  same  per¬ 
centage  change  in  aDF  and  the  region  where  the  opposite  is  true. 
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Pe  WITH  RADAR  SMOOTHING  ONLY  (RS) 

When  only  the  radar  data  are  smoothed,  Eq.  (1)  becomes 


d  = 


1  "  /  -  .2 
R  £(0DFr  ei) 

i  =  i 


:  =  i 


1  "  /  * 

-,-SK 

/=  i 


^DFi  2^2^DFi+^2 


(e2+li) 


2 

n 


X^DFi 

i-  l 


(5) 


This  discriminant  is  identical  to  Eq.  (2).  The  performace  is  therefore  exactly  the  same  as 
when  both  radar  and  DF  data  are  smoothed.  This  equivalence  is  shown  under  a  rather 
restrictive  set  of  assumptions  and  is  not  expected  to  hold  in  a  more  general  case  such  as 
when  the  target  is  in  motion. 


Pe  WITH  DF  SMOOTHING  ONLY  (DS) 

With  DF  smoothing  only,  Eq.  (1)  becomes 


L  <  =  i  J  L  •— 


After  an  unsuccessful  attempt  to  determine  analytically  the  probability  that  this  quantity  is 
less  than  zero,  a  Monte  Carlo  simulation  was  conducted  using  importance  sampling  to  re¬ 
duce  the  number  of  trials  required.  (Importance  sampling  is  briefly  described  in  Appendix 
A.)  To  make  effective  use  of  the  importance-sampling  method,  it  seemed  best  to  use  one 
approach  for  a  <  1  and  another  for  a  >  1.  For  a  <  1  the  parameter  modified  by  importance 
sampling  was  oDF.  This  led  to  an  estimate  of  probability  of  error  of  the  form 
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N,„  n< 

—  Y*  — 

NT  ^aDF 


exp 


U(-d'), 


(7) 


where  0pFfe  is  the  sample  value  of  0DF  on  the  fcth  trial,  NT  is  the  number  of  trials,  o'DF  is 
the  modified  oDF  (the  value  actually  used  in  the  simulation),  and  the  function  U(x)  is  de¬ 
fined  by 


U(x)  =  1,  if  x  >  0, 


=  0,  otherwise. 


The  variable  d'  is  the  computed  discriminant  with  0DF  drawn  from  the  modified  density. 
The  value  of  oR  was  set  to  unity  with  no  loss  of  generality  other  than  the  oR  =  0  case 
treated  separately  below.  Based  on  only  the  assumption  that  Pe  curves  with  DF  smoothing 
would  not  be  drastically  different  from  Fig.  1  (which  turned  out  to  be  correct),  OpF  was 
chosen  to  bring  (nalo’^¥)y/n  to  a  value  (constant  for  each  value  of  oDF  and  chosen  by 
hand)  such  that  Fig.  1  indicated  a  Pe  in  the  range  0.15  to  0.25.  The  results  of  this  simulation 
for  Nt  =  10,000  are  shown  in  Figs.  3a,  3b,  and  3c.  Figures  3d,  3e  and  3f  are  the  results  of  a 
different  simulation;  the  reasons  for  and  description  of  this  simulation  are  as  follows. 

As  o  is  increased  to  values  above  unity,  it  becomes  impractical  to  obtain  an  increase  in 
simulated  errors  by  varying  oDF.  To  see  this,_first  note  that  o  <  1  implies  oDF  >  oR  as 
sketched  in  Fig.  4a.  Most  errors  occur  when  0DF  falls  above  02r  Increasing  oDF  by  a  modest 
amount  will  increase  the  probability  of  this  type  of  error  significantly.  Figure  4b  shows  the 
type  of  situation  implied  by  o  >  1  or  oR  >  aDF.  Here  errors  are  usually  due  to  0Xl  falling 
above  02|.  Increasing  aDF  by  a  modest  factor  will  have  very  little  effect  on  the  probability 
of  error.  Increasing  oD  F  enough  to  significantly  affect  Pe  would  in  fact  lead  to  the  nonzero 
terms  in  Eq.  (7)  (ratio  of  the  true  to  the  modified  probability  density  function  for  0DF) 
being  usually  very  small  but  occasionally  very  large.  This  implies  an  undesirably  large  vari¬ 
ance  for  Pg.  (This  was  in  fact  demonstrated.) 

For  the  special  case  in  which  o  =  aR  =0  the  simulation  is  not  necessary,  because  RDS 
and  DS  become  identical.  (It  is  immaterial  whether  or  not  radar  measurements  are  smoothed 
when  they  are  exact.)  From  examination  of  Eq.  (3)  in  the  limit  as  o  goes  to  0 
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f x .  Vn 


Fig.  3a  —  Performance  with  DF  smoothing  only  (DS):  10,000  trials  with  {^alo^F)^n  set  to  2  and 
a  «  0.25,  with  the  four  curves  being  for  n  «  1,2,  4,  and  8 
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Fig.  3b  —  Performance  with  DF  smoothing  only :  10,000  trials  with  (tijo pF  \l\/n  set  to  2  and  a  -  0.5, 

with  n  •  1,  2,  4,  and  8 
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Fig.  3c  —  Performance  with  DF  smoothing  only:  10,000  trials  with  )\/n  set  to  2.75  and  a  ■  1, 

with  n  ■  1,  2,  4,  and  8 
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Fig.  3d  —  Performance  with  DF  smoothing  only:  10,000  trials  withX'set  to  2.75  and  O-l, 

with  n  “  1,  2,  4,  and  8 
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H  Vn 

Fig  3£  —  Performance  with  DFamoothing  only:  10,000  trials  with  X' set  to  7  and  0-4, 

with  n  -  1,  2,  4,  and  8 
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For  a  >  1  the  quantity  varied  in  the  importance-sampling  formulation  was  p0,  the  true 
separation  in  azimuth  of  the  two  radar  targets.  The  primary  type  of  classification  error  for 
a  >  1  is,  as  described  earlier,  that  for  which  >  02) :•  It  should  be  apparent  from  Fig.  4b 
that  the  frequency  of  this  type  of  error  should  increase  significantly  as  na  becomes  smaller. 
Also,  for  a  <  1  as  in  Fig.  4a  changing  pa  enough  to  significantly  effect  Pe  would  result  in  the 
ratio  of  the  true  to  the  modified  probability  density  function  becoming  usually  very  small 
but  occasionally  very  large.  Therefore,  to  hold  down  the  variance  of  the  estimate  of  Pe 
varying  na  is  restricted  to  a  >  1. 

A  significant  amount  of  computation  was  saved  by  transforming  Eq.  (6)  to  an  equiv¬ 
alent  form 


d'=y(zl~zn+l)  + 


k  =  2 


(8) 


where 


and  where  the  zfe  are  independent  and  Gaussian  with  unity  variance.  All  the  zk  except  z} 
and  zn+ 1  have  zero  mean,  and 


E{z1)  =  Xm^ 


and 


E{zn+ 1  }  =  Xm2  > 


where 


X  =  ny/n= - y/n  , 


DF 


m 


1 


1+7 

207 
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and 


m 


2 


1~7 

207 


The  quantity  X  is  the  quantity  against  which  Pe  was  plotted  in  Figs.  1  and  3.  This  equivalent 
formulation  is  derived  in  Appendix  B.  Varying  pa  in  Eq.  (6)  is  equivalent  to  varying  X 
(which  is  proportional  to  na)  here.  The  computational  advantage  is  obtained  because  only 
two  random  variables  are  affected  now,  whereas  varying  pQ  in  Eq.  (6)  affected  2 n  random 
variables. 

The  error-rate  estimate  is  now 


Pe  Xr^l)2  +(zn*l~  Xm2^ 

T  ft=l 

-  (Z1  -  X'm j)2  -  ( Zn  +  1  -  Ar'm2)2J|  U{-d') , 


where  X'  is  the  modified  value  of  X. 

The  results  of  this  simulation  a re  shown  in  Figs.  3d,  3e,  and  3f.  The  a  =  1  case  was  in¬ 
cluded  to  allow  a  comparison  of  the  two  simulations.  The  smoothness  of  the  curves  in  Figs. 
3d,  3e,  and  3f  does  not  reflect  superior  accuracy  compared  to  Figs.  3a,  3b,  and  3c,  but  in¬ 
stead  reflects  the  fact  that  the  points  on  the  curves  are  not  independent.  Because  the  param¬ 
eter  modified  in  the  importance-sampling  procedure  is  the  parameter  used  as  the  indepen¬ 
dent  variable  in  the  plots,  and  because  the  value  to  which  it  is  modified  (selected  as  before 
for  an  actual  error  rate  of  0.15  to  0.25)  depends  only  on  a,  the  actual  set  of  NT  trials  was 
carried  out  only  once  for  each  combination  of  n  and  a  and  was  not  repeated  for  different 
values  of  X.  The  only  things  that  changed  with  X  were  the  weights  given  to  the  outcomes  of 
the  trials.  In  fact  when  a  calculation  of  d'  at  a  higher  value  of  n  fo’lowed  a  calculation  using 
a  smaller  n  (but  with  other  parameters  unchanged),  only  those  terms  of  d'  corresponding  to 
the  difference  between  the  two  n  values  needed  to  be  calculated.  These  observations  greatly 
reduced  the  amount  of  computation  from  that  required  with  all  points  on  the  curves  in¬ 
dependently  calculated.  They  also  imply  that  both  the  points  within  a  curve  and  the  curves 
within  a  plot  in  Figs.  3d,  3e,  and  3f  are  dependent.  For  example,  comparison  with  the 
curves  for  other  values  of  a  indicates  that  Fig.  3c  may  be  a  better  representation  of  true  Pc 
than  Fig.  3d,  even  though  Fig.  3d  appears  neater. 
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Pe  WITH  NO  SMOOTHING  (NS) 

When  no  smoothing  is  used,  Eq.  (1)  becomes 


d  = 


1=1 


1  n  / 

■in- 

i=  1 


20DFi02  i  +  20DFi0li  +  02  i 


'?<) 


=  n^2  (02j  ~  0li)[(02i  +  0li)  20D  Fi  ]  ' 
i=  1 


With  definition  of  the  variables 


*2  r«i, 


t,,'~  >/2a 


and 


02i  +  0li  20DF i 


W; 


o2  +  2 


d  becomes 


d  = 


2o\/o2  +  2  " 

= - jw 


i=  1 


It  is  easily  shown  that  e(-  and  wi  are  independent  Gaussian  random  variables  with  unit  vari¬ 
ance.  The  calculation  of  their  means  is  postponed  until  needed.  After  the  substitutions 


>l  -  xi 


vi  = 
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and 


Vi 


+  X; 


W:  =' 


d  becomes 


d  = 


Oy/o2  +  2 
n 


(9) 


Solving  for  xi  and  y  ■  in  terms  of  i>(.  and  wi  results  in 


wi  ~  vi 


x-  =- 


s/Z 


and 


w-  +  v- 

i  i 


V* 


It  is  readily  seen  that  the  x  •  and  y  •  are  Gaussian  variables  independent  of  each  other  with 
unit  variance.  The  parenthesized  terms  in  Eq.  (9)  are  therefore  the  independent  noncentral 
chi-squared  [6]  variates  defined  by 


and 


,2,  V-,2 


i=  1 


'2  -  2>? 

i=l 
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with  noncentrality  parameters  a 2  and  6 2  respectively,  where 

«2-£(*w)2 

;=i 


and 


i=l 


Substituting  for  x.  and  y(-  yields 


a 


2 


n/  ~  ^{y,>\  2 

vr— ^ 


and 


where 


and 


£{u>(.}  = 


y/fy'o2*  2 
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Consequently 


y/nn/l  _  1  X 

2  y'  s/c^i) 


The  probability  density  functions  fora:2  andy2  are  [6] 


and 


Rewriting  Eq.  (9)  as 


d  = 


-  x 


2) 


shows  the  probability  that  d  <0  is  the  same  as  the  probability  that  y2  < x2. 


p  x 

Pe  =  j  px2{x2)d(x2)  f  pyt(y2)  d(y2) . 


or 
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Substituting  for  the  density  functions  and  carrying  out  a  change  of  variables  results  in 


x 

v+1e~(x2+a2)l2i^ox)dx J'yv*1e~('y  +b  ^2/„(by)dy, 

0 

where  v  =  (n  -  2)/2  and  n  is  assumed  even.  This  integral  was  evaluated  by  Price  in  Appendix  C 
of  Eef.  7,  which  is  reproduced  in  this  report  as  Appendix  C  (with  comments  and  errata  for 
Appendix  C  being  given  in  Appendix  D).  From  Eqs.  (C-37)  and  (C-50)  the  probability  of 
error  is 


-e(v5'v^)*e 


,-(a2  +  62)/4 


m^m 


m  =0 


(l) 


where 


D 


m 


2~i~v~  i  . 


This  result  is  expressed  in  terms  of  the  Marcum  Q-function,  defined  in  Eq.  (C-4),  and  Im , 
the  modified  Bessel  function  of  the  first  kind.  The  factor  6m0is  the  Kronecker  delta  func¬ 
tion,  which  is  equal  to  one  when  its  subscripts  are  equal  and  is  otherwise  equal  to  zero.  Pe 
can  be  derived  in  closed  form  for  odd  n  by  expressing  it  in  terms  of  the  doubly  noncentral 
F  distribution  [8] ,  but  the  resulting  expression  is  so  complex  as  to  make  its  evaluation  in 
the  present  context  unjustified.  For  the  special  case  n  =  1  the  resulting  Pe  is  identical  with 
that  for  the  other  smoothing  options. 

Pe  was  evaluated  with  the  aid  of  the  techniques  outlined  in  Ref.  9  for  the  evaluation  of 
the  Q-function  and  Ref.  10  for  the  Bessel  functions.  The  results  are  plotted  in  Figs.  5a 
through  5f.  Not  surprisingly,  NS  turns  out  to  be  the  poorest  performer  of  the  four  smooth¬ 
ing  options. 


BOUNDING  PERFORMANCE  AS  RESTRICTIONS  ARE  REMOVED 

Under  the  assumptions  outlined  at  the  beginning  of  this  report  the  best  performance  is 
obtained  (Fig.  1)  when  the  radar  measurements  are  smoothed  prior  to  their  use  in  a  discrim¬ 
inant.  As  long  as  the  radar  measurements  are  smoothed,  performance  is  not  affected  by 
smoothing  the  DF  measurements.  This  is  somewhat  surprising,  because  in  the  absence  of 
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Fig.  5a  -  Performance  with  no  smoothing:  o  =  0,  with  n  =  1,  2,  4,  8,  and  16 


25 


J.  O.  COLEMAN 


/i.v'TT 

Fig.  5f  —  Performance  with  no  smoothing:  a  *  4,  with  n  »  1,2,  4,  8,  and  16 
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smoothing  of  the  radar  data  performance  is  improved  by  smoothing  the  DF  measurements 
(Figs.  3  and  5).  Although  the  probability  of  error  derivations  given  earlier  do  not  generalize 
easily  to  multiple  (more  than  two)  radar  targets,  these  same  conclusions  are  expected  to 
hold. 


With  nonstationary  targets  the  situation  is  expected  to  be  somewhat  different.  The 
simple  averaging  used  for  smoothing  with  stationary  targets  would  have  to  be  replaced  by  a 
significantly  more  complex  procedure.  Examples  of  suitable  procedures  are  a-p  tracking 
[11]  and  Kalman  filtering  [12] .  Suppose  the  radar  and  DF  systems  are  synchronized  in 
time  so  that  when  a  DF  measurement  is  made,  the  radar  system  simultaneously  measures 
the  positions  (in  two  or  three  dimensions  for  best  smoothing  performance)  of  all  of  the 
targets.  Even  under  these  rather  ideal  conditions,  the  quality  of  a  smoothed  estimate  would 
be  reduced  from  the  stationary  target  case,  and  hence  the  remaining  performance  gain  dup 
to  smoothing  would  be  diminished.  The  amount  of  the  remaining  performance  gain  would 
depend  on  the  number  of  degrees  of  freedom  in  the  motion  of  the  targets.  If  the  targets 
were  constrained  to  move  in  a  straight  line  (two  degrees  of  freedom  if  the  motion  is  con¬ 
fined  to  a  plane),  performance  might  be  expected  to  approach  that  of  the  nonstationary 
case  (with  a  large  number  of  measurements).  At  the  other  extreme,  if  no  constraints  what¬ 
ever  were  placed  on  the  motion  of  the  targets,  positions  at  different  instances  in  time  would 
be  completely  independent  and  smoothing  would  contribute  nothing  (fortunately  this  is  not 
realistic).  The  probability-of-error  curves  in  Figs.  1  and  5  therefore  represent  lower  and 
upper  bounds  respectively  on  the  probability  of  error  with  target  motion,  synchronous  mea¬ 
surements,  and  radar  smoothing  only.  With  DF  smoothing  only  and  conditions  otherwise  as 
just  described  the  probability  of  error  would  be  bounded  by  Figs.  3  as  the  lower  bound  and 
Figs.  5  as  the  upper  bound.  (The  assumption  of  independent  unbiased  Baussian  measure¬ 
ment  errors  remains  in  force  throughout  this  discussion.) 

If  the  radar  and  DF  measurements  were  made  asynchronously  and  target  motion  were 
still  permitted  (there  being  no  need  to  synchronize  the  measurements  if  targets  are  guaran¬ 
teed  stationary),  additional  difficulties  would  be  introduced.  Subtracting  an  estimated  azi¬ 
muth  (as  is  done  under  two  of  the  four  options  discussed)  is  meaningful  only  when  the 
estimate  refers  to  the  time  the  measurement  was  made.  Comparing  radar  and  DF  measure¬ 
ments  directly  in  this  context  is  not  possible,  and  some  sort  of  smoothing  or  extrapolating 
procedure  would  be  necessary.  If  only  the  radar  data  were  smoothed,  they  would  have  to  be 
extrapolated  to  the  times  of  the  DF  measurements.  This  type  of  extrapolation  introduces 
additional  error,  as  it  must  be  based  on  inexact  velocity  estimates.  Due  to  this  additional 
error,  the  error  rate  can  never  be  expected  to  equal  the  lower  bound  given  by  Fig.  1.  No 
upper  bound  on  the  probability  of  error  in  this  situation  is  available. 

With  asynchronous  operation  and  target  motion,  smoothing  the  DF  data  without 
smoothing  the  radar  data  introduces  problems.  With  radar  measurements  of  different  targets 
taken  at  different  times,  to  what  time  should  the  DF  data  be  extrapolated?  To  extrapolate 
to  the  time  of  each  radar  measurement  in  turn  could  require  an  unrealistic  amount  of  com¬ 
putation.  Realistic  or  not,  Figs.  3  provide  the  lower  bound  on  error  rate.  No  upper  bound  is 
available. 


Table  2  identifies  the  figures  representing  the  bounds  described  in  this  section. 
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Table  2  —  Summary  of  Figures  Bounding  the  Probability  of  Error 


Target 

Motion 

Allowed? 

Synchronous 

Measurements 

Figure  Numbers  of  Pe  Bounds 

RS 

DS  j 

RDS 

Lower 

Upper 

Lower 

Upper 

Lower 

Upper 

mm 

Irrelevant 

mm 

mm 

3 

3 

mm 

Yes 

WSm 

]9bB| 

3 

5 

BSS 

■si 

No 

mm 

3 

NA 

mm 

CONCLUSIONS 

From  the  discussion  thus  far  it  is  apparent  that  in  the  general  target-motion-with- 
asynchronous-operation  environment,  smoothing  of  the  radar  data  should  be  used.  As  dis¬ 
cussed  earlier  in  this  report,  the  use  of  DF  smoothing  in  addition  to  radar-data  smoothing 
does  not  necessarily  improve  performance.  (Equation  (5)  showed  that  with  stationary  targets 
no  improvement  is  obtained.)  DF  smoothing  can  be  used  in  a  limited  way,  however,  to  reduce 
the  total  computational  burden  in  certain  types  of  environments.  For  example,  the  DF  data 
can  be  preprocessed  by  breaking  up  the  data  into  batches  in  such  a  way  that  all  measure¬ 
ments  within  a  batch  are  taken  over  a  short  enough  time  interval  that  they  can  be  averaged 
together  and  submitted  to  further  processing  as  a  single  measurement.  If  the  environment  is 
such  that  several  independent  DF  measurements  can  be  made  in  a  time  interval  in  which 
total  angular  motion  of  a  target  would  never  approach  the  standard  deviation  of  the  result¬ 
ing  average,  error-rate  performance  should  not  be  degraded.  This  can  be  understood  (at  least 
in  terms  of  the  stationary-target  case)  by  noting  that  if  DF  measurements  are  batched  into 
groups  of  nb  in  the  manner  described,  then  in  Eq.  (4)  oDF  must  be  replaced  by  oDF/\/n^ 

(the  standard  deviation  of  the  average)  and  n  must  be  replaced  by  n/nb.  To  prevent  this 
change  of  n  from  inappropriately  raising  the  variance  of  01  and  0 2  by  nb,  oR  must  be  re¬ 
placed  by  oRly/rTb.  The  effects  of  these  substitutions  exactly  cancel.  The  reduction  of  n  in 
Eq.  (5)  will  more  than  compensate  for  the  computing  time  required  for  the  preprocessing. 

Now  that  the  desirability  of  smoothing  the  radar  data  has  been  established,  how  much 
(or  what  method  of)  smoothing  is  appropriate?  Comparison  of  Eqs.  (1),  (3),  and  (5)  shows 
that  the  performance  difference  between  the  best  options  (RDS  or  RS)  and  each  of  the  other 
is  zero  if  a  =  0  or  n  =  1  and  increases  as  either  o  or  n  becomes  greater.  This  is  consistent  with 
the  intuitively  satisfying  notion  that  the  details  of  the  smoothing  technique  will  not  signifi¬ 
cantly  affect  performance  as  long  as  the  variance  of  the  smoothed  radar  data  is  less  than  the 
variance  of  the  DF  bearing  against  which  it  is  compared.  Both  because  of  radar  data  smooth¬ 
ing  and  because  this  DF  bearing  may  be  the  result  of  averaging  together  a  number  of  DF 
measurements,  the  relationship  of  these  two  variances  may  be  significantly  different  than 
the  relationship  between  the  variances  of  the  measurements  themselves. 
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Appendix  A 

IMPORTANCE  SAMPLING 

Accuracy  in  probability  estimation  by  Monte  Carlo  simulation  is  obtained  by  making 
the  number  of  experimental  trials  sufficiently  high  that  the  event  whose  probability  is  being 
estimated  occurs  many  times.  For  low  probability  events,  the  number  of  trials  required  for  a 
desired  accuracy  can  be  impractically  large.  Importance  sampling  [Al,  A2]  is  a  simulation 
technique  which  can  often  reduce  the  number  of  trials  required  by  modifying  the  probabil¬ 
ity  density  of  the  random  variable(s)  in  the  simulation  to  increase  the  number  of  occur¬ 
rences  of  the  event  of  interest.  The  estimate  is  then  adjusted  in  such  a  way  as  to  remove  the 
bias  this  would  otherwise  induce  in  the  estimate.  The  explanation  which  follows  is  patterned 
after  Mitchell  [Al] . 

First  consider  the  ordinary  Monte  Carlo  simulation  procedure  depicted  in  Fig.  Al.  For 
each  trial  a  Bernoulli  random  variable  y  representing  the  event  of  interest  is  generated  as  a 
function  of  an  input  random  vector  x  characterized  by  a  density  function  p(x).  That  is, 
y  =  F(x).  The  variable  y  is  equal  to  one  if  the  event  of  interest  occurred  and  zero  if  it  did 
not.  The  probability  estimate  is  obtained  as  the  average  of  y  over  many  trials. 

For  the  simulation  using  importance  sampling  shown  in  Fig.  A2  the  density  p(x)  is  re¬ 
placed  with  a  modified  density  pm  (x)  which  will  cause  the  event  of  interest  to  occur  more 


PROBABILITY 

ESTIMATE 


Fig.  Al  —  Ordinary  Monte  Carlo  simulation 


Fig.  A2  —  Simulation  with  importance  sampling 
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frequently.  The  processor  output  is  multiplied  by  a  weighting  function  w(x)  to  compensate 
for  the  change  from  p(x)  to  pm  ( x ).  The  weighting  function  must  be  chosen  so  that 
E{y}  =  E{w(x)ym}or 


f F(x)p(x)  *-/' 


F(x)w(x)pm(x)  dx 


Therefore  the  weighting  function  must  be 


w(x)  = 


P(x) 
Pm(x)  ' 


The  weight  to  be  applied  for  a  particular  trial  is  just  the  ratio  of  the  original  to  the  modified 
density  function  evaluated  at  the  specific  value  of  x  used  for  that  trial. 

The  variance  of  the  estimate  can  be  expressed  in  terms  of  the  first  and  second  moments 
of  the  estimate.  Because  the  estimate  has  been  designed  to  be  unbiased,  the  first  moment  is 
equal  to  the  probability  being  estimated.  Obviously  in  any  case  of  real  interest  this  will  not 
be  available  in  a  suitable  mathematical  form  or  there  would  be  no  need  to  perform  the  simu¬ 
lation.  The  second  moment  tends  to  be  at  least  as  elusive.  In  the  absence  of  analysis  of  the 
variance  of  the  estimate,  the  choice  of  a  modified  density  pm  (x)  becomes  heuristic.  It  can 
be  shown,  however,  that  in  cases  involving  simple  (analytically  tractable)  functions  p(x)  and 
F(x),  importance  sampling  can  dramatically  improve  estimation  accuracy  [Al,  A2] . 


REFERENCES 

Al.  R.  L.  Mitchell,  “Importance  Sampling  Applied  to  Simulation  of  False  Alarm  Statis¬ 
tics”,  to  be  published  in  IEEE  Trans,  on  Aerospace  and  Electronic  Systems  (around 
January  1981). 

A2.  V.  G.  Hansen,  “Importance  Sampling  in  Computer  Simulation  of  Signal  Processors”, 
Comput.  and  Elect.  Engrg.  1  (No.  4),  545-550  (1974). 


35 


Appendix  B 

DERIVATION  OF  EQ.  (8) 


Equation  (6)  can  be  rewritten  using  vectors  rather  than  summations  as 
d  =  “l|0DF_  1|  '  “||0DF~  01 1| 

by  defining  the  n-vectors 


o 

9  21 

o 

II 

^DF 

a 

»-H 

*12 

*  ®2  " 

9  22 

1 

a 

_V 

9  2  n 

The  vectors  0DF,  ®i>  and  02  are  independent  Gaussian-distributed  random  vectors  with 
covariances  M,  o2I,  and  a2/  respectively,  where  M  is  an  n-by-n  matrix  with  every  element 
equal  to  1/n  ancf  /  is  the  n-by-n  identity  matrix. 

By  defining  (for  this  appendix  only)  a  partitioned  vector 


(®DF_  0_2) 
(®DF  ~®l) 


and  a  partitioned  matrix 


Q  = 


the  discriminant  equation  can  be  written  as 

\ 


nd  =xtQx  . 
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The  Gaussian  random  vector  x  has  a  mean 


E\x\ 


6 


where  c  is  defined  to  be  an  n-vector  whose  elements  are  all  unity.  The  covariance  of  x, 
denoted  by  K,  is  found  to  be 


a2/  +  M 

M 

M 

a2I  +  M 

From  this  point  on,  linear  transformations  on  x  are  used  to  proceed  toward  the  desired  dis¬ 
criminant  form. 

It  is  simple  to  demonstrate  that  the  matrix  K  is  positive  definite  as  long  as  o  is  nonzero. 
Because  it  is  both  positive  definite  and  symmetric,  it  can  be  factored  into  symmetric  posi¬ 
tive-definite  square  roots: 


K  =  SS. 


It  is  straightforward  to  verify  that  S  is  given  by 


I 


I 


Because  S  is  positive  definite,  its  inverse  must  exist  and  is  obtained  as 

1  /I  1 


■  1  _ 


a__-2___v/°__+2_  _ 
2 \°  +  2  / 


2  \a  ,/a2  +  2 
__  -----  f--,-- 

-I-  -  -:  )M 

0  2\° 
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and  used  to  transform  x  to  a  new  variable 

y  =  S-1x  . 


The  covariance  of  y  is 

£j(y- E{y})(y-£{y})Tl  =  S-1Ej(x- E{x})(x-E{x})T[s-1 

=  S~1KS~ 1 
=  / 

and  the  mean  of  y  is 

£M 


The  discriminant  equation  can  now  be  written  in  terms  of  the  transformed  variable  y  as 

nd  =  xtQx  =  yTSQSy  =  yTTy  , 


where 


T  =  SQS . 

Performing  the  multiplications  and  simplifying  yields 


al  +  (  y/o2  +  2  -  a  ) 

Im  i 

i 

0 

0 

\-aI-  | 
1 

[•Jo2  +  2-  o  j 

|  Af 
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The  next  step  in  the  derivation  requires  that  T  first  be  diagonalized.  Matrices  U  and  D 
are  required  such  that 


T  =  UtDU  , 


where  U  is  orthogonal  (UT  =  U~  1 )  and  D  is  diagonal  (all  off-diagonal  elements  are  zero). 
The  matrix  U  satisfying  these  requirements  is  in  this  case  of  the  form 


U  = 


where  F  is  any  n-by-n  orthogonal  matrix  whose  first  row  is  composed  of  identical  elements. 
It  is  not  necessary  to  specify  F  further  in  order  to  complete  the  derivation.  For  convenience 
the  value  in  the  first  row  will  be  assumed  positive,  which  implies  a  value  of  1  jy/n.  A  familiar 
example  of  such  an  F  matrix  is  a  discrete-Fourier-transform  (DFT)  matrix. 

Solving  for  the  diagonal  matrix  D  results  in 


v/l  +  2  la2  0 
0  1 

0  0 

0  0 

1 

0  •••  0  j 

0  -  0  j 

1 :  o  i 

...  o  1  ! 

0 

2/a2  0  0-0 

i  0 

-1  0  ■  0 

0 

o 

0  -1  . 

.  0 

1  0 

1 

0  0  -1 

A  new  random  vector  can  now  be  defined  as 


2  =  Uy  , 


and  the  discriminant  equation  becomes 


nd  =  yTTy  =  (UTz)T(UTDU)(UTz)  =  zTDz  . 
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The  covariance  of  the  new  variable  is  the  same  as  the  old : 

cov  {  z }  =  cov  {Uy}  =  U  cov  {  y)  UT  =  UIUT  =  I . 

The  mean  of  z  is 

E{z}  =  UE{y) 


If  z  is  expanded  to 


the  discriminant  equation  can  be  written  in  nonmatrix  form  as 
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where  the  zk  are  independent  Gaussian  random  variables  with  unit  variance.  Defining 


d'  = 


nd 


a 


2 


and  noting  that  the  means  of  the  z  values  here  match  those  in  the  text  completes  the 
derivation. 


3 
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APPENDIX  C 


APPENDIX  C  [Reproduction  of  Appendix  C  of  Ref.  7] 
REDUCTION  OF  A  DOUBLE  INTEGRAL  INVOLVING  BESSEL  FUNCTIONS* 

R.  Price,  MIT  Lincoln  Laboratory 


I.  INTRODUCTION 

Maximon2<*  has  studied  the  following  integral,  which  for  v  =  0  arises  in  trying  to  find  the 
probability  that  the  envelope  y  of  a  sine-wave-plus-narrow-band-noise  process  is  less  than  r 
times  the  envelope  x  of  another  such  process,  the  noises  being  assumed  independent  and  both  of 
unit  variance. 

P  (a,b,r)  =  a^b"1*  C  x"^+1  exp  [—  (x2  +  a2)/2)  l^fax)  dx 
Jo 

•  C”  /+1  e*pMy2  +  bV2)  Vby>  dy  <c-u 

Jo 


Here  lj.z)  is  the  modified  Bessel  function  of  nonnegative  integral  order  v,  and  a  and  b  are 
parameters  which,  for  v  =  0,  are  respectively  equal  to  the  amplitudes  of  the  sine  waves  in  the 
processes  whose  envelopes  are  x  and  y.  (See  Rice3^  for  the  derivation  of  the  envelope  proba¬ 
bility  density  functions  contained  in  Eq.  (C-l),  with  v  -  0.)  Maximon  obtains  P  in  terms  of  a 
Neumann  series  of  modified  Bessel  functions: 


P_„(a,b,r)  =  (^)Z>'exp(-y  (a2r2  +  b2)/(l  +  r2)]  £ 

m=0 


,  ar.Hm  , 

m  '  b  '  V+m 


/  abr  \ 

'777*' 


(C-2) 


where 

r2/(l  +  r2)  (if  m  =  0) 

c  = 

m  1  (if  m  >  0)  ( v  an  integer  >  0)  (C-3) 

The  series  (C-2)  can  be  considered  in  closed  form  in  terms  of  Lommel’s  functions  of  two  vari¬ 
ables,31  but  tabulations  by  which  the  series  may  be  evaluated  exist  chiefly  in  terms  of  the 
"Q-function"  examined  by  Marcum.32  This  function  is  the  cumulative  probability  distribution  of 
the  envelope  of  a  sine-wave-plus-narrow-band-noise  process: 

Q(u,v)  =  C  x  exp(-(x2  +  u2)/ 2]  I  (ux)  dx 

Jv 


=  1  -  exp(-(u2  +  v2)/2]  Yj  <-7>m  (uv) 
m=l 


=  exp[-(u2  +  v2)/2J  £  (_7)m  !m  (uv> 

m=0 


(C-4) 


*  Jon**33  hoi  recently  and  independently  reduced  the  double  integral  [Eq.  (C-6)  J  for  v  -  t  =  r,  obtaining  the  par¬ 
ticular  com  of  the  solution  [Eq.(C-35)|  where  m  =  v  =  1  [given  alio  by  Eq.(C-41)). 
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The  first  series  is  due  to  Bennett,  and  the  second,  which  is  connected  with  the  reciprocation 
properties  of  Loramel's  functions  of  two  variables,31  is  quoted  by  Helstrom  23 

In  this  Appendix  we  deal  with  a  double  integral  similar  to  Eq.  (C-l),  which  arises  in  finding 
the  probability  that  the  sum  Y  of  a  number  (v  +  1)  of  Squared  envelopes  of  sint -wave-plus- 

narrow-band-noise  processes  is  less  than  r2  times  another  such  sum  X  where  all  the  noises  are 

*  14 

independent  and  have  unit  variance.  The  probability  density  function  of  X  is 

p(X)  =-i-(X/o)l'//2  exp  [—  (X  +  a)/2]  1^  */xZT  (C-5) 

and  that  of  Y  is  identical  to  Eq.  (C-5)  with  0  in  place  of  a.  Here  a  and  0  are  proportional  to 
the  sums  of  the  sine-wave  powers  (squared  amplitudes)  in  the  processes  yielding  X  and  Y,  re¬ 
spectively.  The  double  integral  for  the  probability  that  Y  is  less  than  r2  times  X  is  thus,  after 
making  the  changes  of  variables  X  =  x2,  Y  =  y2  and  a  =  a2,  0  -  b2: 

/~»oo 

P^a.b.r)  =  (ab)"1*  \  x^+1  exp[-(x2  +  a  )/2)  1^  (ax)  dx 

•  C  y*'+1  exp  [—  (y2  +  b2)/2]  1  (by)  dy  (C-6) 

Jo 


H.  DEVELOPMENT  OF  A  RECURSION  RELATION  FOR  P„(a,b,  r) 

The  first  step  in  attacking  Eq.  (C-6)  is  to  integrate  by  parts: 

exp[(a2  +  b2)/2)  (ab)*' P^a.b.r)  =  -  exp  [—  x2/ 2]  (x^Ux)  ^  (  )  dy] 


rx  /2  -i- 

dx 


K  *.*•«>  f 


,(by)  dy 


dx 


(C-7) 


where,  for  v  >  1,  the  first  term  on  the  right  of  Eq.  (C-7)  vanishes.  We  now  make  use  of  the 
34 

Bessel-function  property 


which  yields 


1  d  ,  v,  ,  ,,  v- 

z  dl  Iz  Vz)l  *  z 


,(z) 


(C-8) 


^  (x^fcx)]  =  cx^^cx)  (C-9) 

Expanding  the  derivative  in  Eq.  (C-7)  and  applying  Eq.  (C-9),  we  obtain 

exp((a2  +  b2)/2](ab)1' P^(a,b,  r)  =a  y  xv  exp [- xZ/2] I^^ax)  j ^  y*'+1  exp(-y2/2) I^fby)  dy]  dx 

+  r,'+2y  x2”*1  exp|-x2(l  +  rZ)/2]yax)ybrx)  dx  .  (C-10) 

Again  integrating  by  parts  in  the  inner  integral  of  Eq.  (C-10),  and  using  Eq.  (C-9), 


*  In  farms  of  classical  statistics,  this  is  equivalent  to  finding  the  probability  distribution  function  of  the  ratio  of 
the  lengths  of  two  (2v  +  2)-dimensional  random  vectors,  all  of  whose  components  are  mutually  independent  and 
Gaussian  with  unit  variance,  and  whose  two  mean  vectors  are  of  specified  lengths.  It  appears  that  the  method 
of  solution  followed  here  could  be  generalized  to  vectors  of  different  dimensionalities. 
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C  yl<+1  exp[-y2/2]  I  (by)  dy  =  -expf-y2/2)  yl'l  (by)  I 

+  b  y  yv  exp[-y2/2]  I[/_1(by)  dy 


(C-ll) 


which,  when  substituted  into  Eq.  (C-10)  yields,  for  v  1, 

exp((a2  +  b2)/2](ab),'P  (a,b,r)  =  ab  C  xv  exp[-x2/2)I  ,(ax)  dx  f  yv  exp[-y2/2)I  j(by)  dy 

Jo  Jo 

+  rv+Z  C  x2l/+1  exp  [-  x2(i  +  r2)/ 2 ]  I  (ax)  I  (brx)  dx 

-ar"  C  xZv  exp[-x2{i  +  r2)/2]Ip_  (ax)  I  (brx)  dx  .  (C-12) 

In  order  to  cast  Eq.  (C-12)  in  a  more  symmetric  form  for  later  convenience,  we  note  that,  inte¬ 
grating  by  parts  and  using  Eq.  (C-9)  once  again,  with  v  >  1, 

„+2r  2v+i  -  '  -i,+z--zv 


C  xZv+l  exp(-x2(l  +  r2)/2]  I„(ax)  I  (brx)  dx  =  -  - - exp[~x2(l  +  r2)/2]  1  (ax)  I  (brx) 

J0  v  v  1  +  r 

+  -—j  j"”  expf-x2(l  +r2)/2)  {[x\(ax))  (x^ybrx)]}  dx 

v+Z  p"°  ,  ,  _ 

=  \  x21’  exp  (-x2(l  +  rZ)/2]I  ,(ax)  L(brx)  dx 

1  +  r  Jo 

+  -^4C  xZv  exp(-x2(l  +  r2)/2)I  (ax)  I  ,(brx)  dx  . 

1  +  r  Jo 


(C  —  1 3 ) 


Combining  Eqs.  (C-12)  and  (C-13),  we  have 

P^(a,  b,  r)  =  (ab)1"I/ y  xV  expf-xVzjI^lax)  dx  J  yv  exp[-y2/2]  IJ/_1(by)  dy 


o 

.1  .  ^ 


+  ~ — *  k  |r,/+3  3i“vb*~v  C  x2^  exp[-x2(l  +  r2)/2]l  (ax)  I  .(brx)  dx 

1  +  r  *  v 

-rl/a1"I/b"v  C  x2v  exp[-x2(l  +  r2)/2] I^fax)  I^(brx)  dx} 

^  n  * 


thus  giving  the  recursion 

P  (a,b,r)  =  P  (a.b.r)  +  exPl~(^ 

1  1  +  r 

•  |rv+:3  a'^b1"*1  y  x2v  exp[-x2(l  +  r2)/2]l(,(ax)  lj/_1(brx)  dx 

-rva.i~vb~v  y  x2v  exp(~x2(l  +  rZ)/2JIl,_1(ax)  I^fbrx)  dx)  .  (C - 14) 
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Equation  (C-14)  is  a  reduction  formula  by  which  we  can  obtain  the  general  result  if  we  know 

29 

P0(a,b,r)  as  given  by  Maximon.  It  is  necessary,  however,  to  evaluate  the  pair  of  integrals 
in  Eq.  (C-14),  and  this  is  where  the  complication  lies*  The  integrals  are  solved  by  using  another 
Bessel-function  property35  that  is  similar  to  Eq.  (C -8> : 


,  i  d  -m  ,  -v-m  .  ,  . 


f  C —15) 


which  leads  to 


<7  fa*  (cx)'m  lm^  -  xl"m  c""m  Iv+m<«> 


and  as  a  special  case 


l0{cx)  =  xvc-v  I„(cx) 


Therefore, 


(C— 16) 


(C-t?) 


a-W'*' 


y  xZv  exp[-x2(i  +  r2)/2]Il,(ax)I(,_1(brx)  dx  = 


br  d(br) 


v-1 


xj  x  exp[-x2(l  +  r2)/2]IQ(ax)  IQ(brx)  dx 


(C-18) 


and  similarly  for  the  other  integral  in  Eq.  (C-14),  a  and  (br)  simply  being  interchanged. 

.  * 


The  integral  in  the  right  member  of  Eq.  (C-18)  is 

C  x  exp  (— x2(i  +  ra)/2]  I  (ax)  I  (brx)  dx  =  (1  +  r2)"1  exp  [--- --b-^rZ  1  I  1-5^1  .  (C-19) 

Jo  o  o  1  2(1  +  r2)  1  °  1  1  +  r  ‘ 


Combining  Eqs.(C-14),  (C-18)  and  (C-19), 

.2  .  ,_Z 


l/z'  4^)-’  afei'-1 


•  ~p  iy^-1  ■„  1-^1} 

'  ^ - *■ ■„ Iffy)  *  I^D  «-»> 


so  long  as  v  >  1. 


III.  FINAL  SOLUTION  FOR  P„(a,b,  r) 

Iterating  the  recursion  formula  (Eq.  (C-20)]  and  reverting  to  the  variables  a  =  a2,  (3  =  b2, 
we  obtain 


*Such  infegrali  are  treated  by  Watjon  (Ref.  31 ,  pp.  395-396)  in  connection  with  Weber's  second  exponential  in¬ 
tegral,  of  which  Eq.  (C-19)  is  an  example;  with  regard  to  integrals  of  the  type  appearing  in  Eq.  (C-14)  Watson 
seem*  unduly  pessimistic. 
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P  <^.  n (p.  r)  =  PjVZT.  VJ,  r)  +  «xpJ-  <“  *  m 
v  °  4<r  4  r'V 


j=i 


•  |(r  —  r’1)  Io(Va^r/(l  +  r2)]  4  [r2^  "r*2/fj  ^[VSJr/d  4  r2)  ]|  (C-21) 

which,  expanding  according  to  the  binomial  formula  for  the  derivative  of  a  product,  differentiating 
exp[(a  4  0r2)/2(l  4  r2)),  and  returning  to  the  variables  a,  b,  yields 


P  (a,  b,  r)  =  P  (a,  b,  r)  4 


exp  I—  -- ---  g-  l  v  j-1  j-1 
1  2(1  4  r  )  *  v  V  V 


_ -2<V3r  £  £  £  (1  4  r2)k+i+2-2;  r2j-2/-2 

(r  4  r  ) 


where 


j=l  k=0  1=0 

■  ft1)  f,')  4  {«—-*>'„  |^| 

(D-.il.-iJI  (if  m  <  0  or  m  >  n.  (")■») 


(C-22) 


(C-23) 


is  the  binomial  coefficient. 

We  now  need  to  find  the  result  of  operating  on  the  term  in  the  braces  of  Eq.  (C-22).  Using 
Eq.  (C-16)  with  c  =  b,  x  =  ar/(l  4  r2),  v  =  t  and  m  =  0  or  -1  [note:  Im(z)  =  I_m(z) ),  leads  to 

l-rr^l)  ■(777)'^'-' I777I 


■i  4>''o  (7771) '  (777)' 


(C-24) 


and  further  applicat ion,  this  time  with  c  =  a,  x  =  br/(l  4  r  ),  v  =  k,  and  m  =  1  -  l  or  -t  yields 

,k-f4i  , 


4 &>k  <i  4*'  {t  >,  I777D  •  (777)  k"  W,  [777I 
4  4'"  <i£>‘  {■. I777I}  ■  (777) k"  I777I 

Similarly,  we  find 

,1  d  ,k,l  d  ,1  fa  ,  [  abr  1\  _  /  r  \  k4l  .  ./-k41  I  abr  | 

(a  a?1  <b  db>  lb  h  I777J/  -  (777)  <a/b)  li-M  I777I 


(C-25) 


(C-26) 


(C -27) 


Substituting  Eqs.  (C-25)  through  (C-27)  in  Eq.  (C-22),  and  interchanging  the  dummy  variables 
k  and  t  in  the  part  of  Eq.  (C-22)  corresponding  to  the  term  r'2  (a/b)  I j  [abr/(l  4  r2)] ,  we  are 
led  to 
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exp  -  — - 2— I  v  j-1  j-l 

P„<a,b.r)  =  P  (a.b,r>  + - ^  -  %  (r  +  r‘ V2j  £  £  (^‘l  ft1) 

<r  +  r  1  j=l  k=0  1=0 

.  {(r  _  r-*,  (£)“-'  ik  |  (-^1  +  [r(  >^)k-'+1  -  r-1  (^)k-,+  1  1 


77^1) 


(C-28) 


Letting  m  =  k  —  I  +  1,  and  noting  that  k  and  hence  m  can  range  from  —  »  to  +  “  without  changing 
the  value  of  Eq.  (C-28)  because  of  the  behavior  of  ft1)  as  defined  by  Eq.  (C-23), 

.  2  2  A .2  , 

(  a  r  +  b  t 
exp  - - -j—  +•> 

Pw(..b.r)-P0(a.b.r)+ - ^  ^  1  £  {[r(^)m  -  r*1  (fp)"1] 


+  <r  -  r‘1» l 777)}  l  <r  +  r‘1>'2j  I  <j'/>  LV/.i)  •  <c-29> 

1  +  r  j=i  1=0 

We  now  inquire  about  the  final  sum  in  Eq.  (C-29).  We  have 

£  (2j'2)  xq  =  (1  +  x)2j-2 


=  (1  +  x)^'1  (1  +  x^*1 

■  s'  ‘i  • 

1=0  n=0 

Comparing  powers  of  x,  with  n  =  m  +  I  —  1,  we  deduce 

jfi'iujj-isa 

1=0 

By  Eq.  (C-31)  and  the  behavior  of  ( as  defined  by  Eq.  (C-23),  Eq.  (C-29)  now  becomes 


(C-30) 


(C-31) 


.2.2  .  v2. 


Pv(a,b,  r)  =  PQ 


exP(-^^4_|  v  5 

o<a’b-r)+ - l  (r+r-V2J  £  (2^) 

(r  r  1  j=t  m=-j+2 

■„  I777I  *  i- ' •'*  V.,  I777I} 


-4 i£?fl  ifahhizl 


(r  +  r  ) 


(r  +  r  ) 
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TT^hrW* 

j=2  'm--j+2 

+  r  n  =  -j+l  ’  r  '\ 


exp  I-TTTVI 

P,(a.b.r)-P0(..b.r)  + - —  yEo - 

(r  +  r  )  (r  +  r  ) 


(r  -  r'1) 


-^Al0[hl+  l  Ir  +  r-V^ffltr-r^iyhl 

(r  +  r  >  j=2 


.  br  >  m 
'  a  ' 


*  Z  Z  (H5J.3--*(^)I 

j=2  m  =  l 

+  (r  -  r”1)  £  (r  +  r'V^gflyh] 

j=2 


Mr-r'1)  ^  <r  +  r-V2j  £  (ffl  K^*" *  <&>"]  ^fh] 
j=2  n=i 

where  we  have  used  the  fact  that  I  (z)  =  I  (z),  and  h  =  abr/(l  +  r2).  Gathering  terms, 

m  -m 


P  (a,  b,  r)  =  P  (a,  b,  r)  + 


.  2  2  ,  ,2 , 

exp[--^-^ - j-l  .  v  /2j-t\ 

-  --2(4.^  (r-r'Syh)  £  -Li* 

(r  +  r  *)  |  °  j  =  1  (r  +  r  VJ 


+  (r  +  r-V2f^---5z)  IJh] 

br 

*h-r-'l  Z  lr  ^-V21  l  (/if,)  *  (&>"l  ■„!»] 

j=2  n=l 


Z  (r.r'V2'  Z 

j=2  m  =  l 

"  (r_1  (2ji^)  -  r(?+^)l  (¥F,m)  'm^l 
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I  a2r2  +  b2l 
exp - 5— 

=  Po(a.b.r)  + - |(f  ~  r'1)  Ip[h)  £ 


j=l 


(V) 

(r  +  r”1)2^ 


j=l 


S  V=i  £  (|f^)-r-*(^)|<^V 

m=l 


(C  —33) 


Rearranging  the  double  summation  in  Eq.  (C-33)  so  that  summing  on  m  occurs  first,  we  have, 
finally, 

v 


P„(a,b 


i  (<*>"  i 

1=1  +  r  '  j=m 


x  l-'W -«)!'■■ 

j=m  n 


l  (  '™m°)  *  r'V2i-> 


j=m+<5 


v  v  4  2  ‘ 

V  „r/br,m  ,a,m,  T  |  abr  1  v  (2j  -  1)!  (r  +  r  )  J 

h  -(b?)  1  rm  lTT~Z|  L  “V+"^?ir=mTT— 

m=l  1  +  r  j=m 


(C-34) 


Equation  (C-34)  is  valid  for  v  ^  1 ,  a  ^  0,  and  b  ^  0,  and  6  is  the  Kronecker  delta -function, 

mo 

6 _ -  t  for  m  =  0,  6 _ =  0  otherwise.  When  r  =  1,  Eq.  (C-34)  reduces  to 

mo  mo  *  * 

P^a.b,  l)=Po(a,b,  1)  +  exp  (-(a2  +  b2)/4] 


m=i 


(C-35) 


]=m 


There  is  an  alternate  expression  for  the  final  summation  in  Eq.  (C-35)  that  may  sometimes 
be  preferable.  Using  the  identity 


V  (2j  -  1)1  2'Zj  _  (2m)*‘  V  (v  4  j)!  l~v~i 

“  n+Wnj  -  m).’  (v  +  m).1  ^  (j  -  m)!  *m 

j=m  ]=m 


(C-36) 
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we  obtain  in  lieu  of  Eq.  (C-35) 

Pyfa,  b,  1)  =  P0(a,b,  1)  +  2'1'"1  exp[-(a2  +  b2)/4) 


£  UTT^y.  f(T'm-<¥)mi1m(ab/2)  2 


{V  ±  ■))■'  Z~ 
(j-m): 


(C-3  7) 


To  prove  Eq.  (C-36),  it  suffices  to  establish  that  it  is  true  for  v  =  m,  which  can  be  seen  im¬ 
mediately  because  the  summations  of  Eq.  (C-3 fc)  are  then  just  single  terms,  and  that  the  difference 
in  the  right  side  between  the  summation  for  v  =  M  +  1  and  that  for  m  =  M  is  equal  to  a  like  differ¬ 
ence  in  the  left  side.  Taking  these  differences,  we  have  to  prove 

(2M  4-  1)1  2~2M'2  (2M  +  2)!  2~2M~3 

(M  +  1  +  m).'  (M  +  1  -  m).'  ~  m(M  +  1  +  m).'  (M  +  1  -  m).' 

M  M 

2-m-^  (M  +  i).'2~J(.i  -  M  -  2m  -  1)  ,r 
m  L*  (M  +  1  +  m)!  (j  —  m)! 


which  becomes,  multiplying  Eq.  (C-38)  through  by  2." 
between  the  single  terms 


m(M  +  1  +  m)!  and  taking  the  difference 


1V1 

V  (M  +  ,)!  2'J(M  +  1  +  2m  -  i)  (2M  +  1) !  2" 
^  =  ~  (M  -  mTT 


The  left  member  of  Eq.  (C -39)  may  be  written 


l  Vi 

M  +  J+1)-2(j-m,H  l 


-  I  ij^r^  (c-40) 

j=m+i 

and  changing  the  subscript  j  to  j'  =  j  —  1  in  the  second  summation  of  the  right-hand  side  of 

Eq.  (C-40),  we  obtain  a  sum  identical  to  the  first  sum  on  the  right-hand  side,  except  that  it  reaches 

only  to  M  -  1  rather  than  M.  Thus  Eq.  (C-40)  is  equal  to  (M  +  j  +  1).'  2"ty(j  -  m)!  for  j  -  M,  and 

Eq.  (C -39)  is  thereby  proved,  substantiating  Eq.  (C-36). 

Although  Eqs.  (C-35)  and  (C-37)  are  general  formulas,  it  may  on  occasion  be  preferable  to 

express  the  solution  in  terms  of  P  (a,b,  1)  and  a  single  pair  of  Bessel  functions,  by  reducing 

o  36 

the  sum  of  Bessel  functions  through  Bessel-function  recursion  relations.  We  have  done 
this  for  v  of  1  through  5,  inclusive,  and  find  that  the  minimum  number  of  algebraic  terms 
is  obtained  when  reduction  is  made  in  terms  of  I^fab/2)  and  1^  j(ab/2).  However,  the  number 
of  algebraic  terms  appears  to  grow  like  t  +  [v(v  —  i)]/2,  so  that  it  is  difficult  to  go  much  beyond 
v  =  5: 

Pjfa.b,  1)  =  Po(a,b,  1)  +  2'4(b2  -  a2)  exp[-(a2  +  b2)/4)  q'1  Ijlq)  (C-41) 


50 


NRL  REPORT  8443 


P2(a,b,l)  =  Po(a,b,  1)  +  2*7(b2-a2)  exp[-(a2  +  b2)/4) 

•  [i2q_1  ij(q)  +  q'2  *2(q)  (s) ]  (C-42) 

P3(a.b,  1)  =  PQ(a,b,  1)  +  2*10(b2  -  a2)  exp|-(a2  +  b2)/4J 

•  [16q-2  I2(q)  (29  +  a)  +  q'3  I3(q)  (28p  +  s2)]  (C-43) 

P4(a,b,  1)  =  P0(a,b.  1)  +  2'13  (b2  -  a2)  exp  [-(a2  +  b2)/4] 

•  [4q‘3  I3(q)  (195  ■  25  +  69  •  22S  +  60p  +5s2) 

+  q"4  I4(q)  (65  ■  24p  +  44ps  +  s3 ) ]  (C-44) 

P5(a,b,  1)  =  PQ(a,b,  1)  +  2'16(b2-a2)  exp[-(a2  +  b2)./4) 

■  [8q'4  l4(q)  (843  ■  28  +  8?  •  2?s  +  3104p  +  268S2  +  52ps  +  3S3) 

+  q"5  I5(q)  (843  ■  26p  +  87  •  25ps  +  64pS2  +  496p2  +  s4)]  (C-45) 

2  2  2  2 
where  p  =  a  b  ,  q  =  ab/2,  and  s  =  a  +  b  . 

TV.  SIMPLE  DERIVATION  OF  PQ(a,  b,  r) 

For  completeness  we  present  here  a  short  derivation*  of  the  relation  between  Po(a,b,r) 
and  Marcum's  Q-function,  which  is  quoted  in  Eqs.  (36)  and  (37)  and  which  can  be  obtained  from 
Maximon's  result  (Eqs.(C-2)  and  (C- 3 ) J  and  the  second  series  of  Eq.(C-4)  of  this  appendix.  The 
following  derivation  given  here  appears  to  be  considerably  less  complicated  than  that  employed 
by  Maximon. 

Referring  to  Eq.  (C-6)  of  this  appendix  and  to  Eqs.  (28)  and  (32)  of  the  text  (generalized  from 
2  2 

r  =  1),  with  a  -  a  ,  p  -  b  , 

.  p"°  _  pXr2 

PQ(  •Ja.  r)  =  \  \  exp(-(X  +  a)/2]I  (VXa)  dX  \  exp[-(Y  + /3)/2)  I  (nTy^)  dY  (C-46) 
The  Laplace  transform  of  PQ('/cr ,  *J]f,  r)  with  respect  to  is 

L(a,r,  s)  =  C  P  ('/a, '/if,  r)  e~^s  d/J 
Jo  0 

j  p"  _  p  Xr2  p°°  _ 

=  -4!  exp[-(X  +  a)/2]I  WXa  dX  \  exp(-Y/2)  dY  \  exp[-0(s  +  y]  I  (\fY0)  d0 

Jo  Jo  Jo  ‘  0 

=  4-  C  exp[-(X+c«)/2]l  (n/X^)  dX  CXr  exp[-Y/2  -t  Y/{4e  4  2)]  dy 

Jo  0  Jo  (s  +  j) 

*  The  Laplace  transform  method  given  here  was  apparently  anticipated  by  J.E.  Storer  in  unpublished  work  performed 
In  1960.  (See  the  remarks  of  Jones33  which  introduce  his  Appendix  B.) 


J.  O.  COLEMAN 


J  e-(X+a)/2  |,  _  expf_sxr2/(l  +  2s) )J  dX 


exp  (  - 


L(a,  r,  s)  =  —  - 


{or  I  (  or  1 

- rr  expi — i - 1 - ~i - z t 

2(1  +  r  ) J  .  1  4(r  1  1)  [s  +  2  (1  +  r  )  ]  J 


(1  +  r2) 


s  +  2'1  (1  +  r2)'1 


(1  1  r  )  ']■ 


exp- 


f  or  1  f  err2  1 

1  - rr  expi - 2 - 2 - 7i - TTi  ( 

L  2(1  +  r  ) J  l4(r  +  lpls  +  2  1  (1  +  r  )  1 1  •* 


2(1  +  r  ) 

where  we  have  used  the  Laplace  transform  formula^7 


4(r  4-  1)  [s  +  2  (1  +  r  I'1)' 

s[s  +  2*1  (1  +  r2)'1) 


r°o 

\  exp(-yU)  I  (2x/r)U)  dU  =  y'1  ex p  ( tj /y) 


(C-47) 


(C-48) 


Examining  Eq.(C-47),  we  see  that  it  is  the  Laplace  transform  of  a  sum  of  three  terms;  the 
first  is  unity,  the  second  is,  by  Eq.  (C-48), 

2 


exp 


-/  -"-r _ l 

-llil  +  A-’  .  cxp  / - ^  \  i 

(1  +  r  )  *•  2(1  +  r  ) J  ° 


+  r2)] 


and  the  third  is 


exp  [ - — — y  1  a 

- -  -2ii4  r-LJ  •  f  exp  { - 0-y  }  I0  K5flr/(1 

2(1  +  r  )  Jo  1  2(1  +  r  ) J 


+  r  ) )  d/7 


(C  -49) 


If  the  upper  limit  of  the  integral  in  Eq.  (C -4 9)  were  infinity,  the  entire  expression  would  be  unity, 
again  by  Eq.  (C-48),  so  that  we  conclude 


P  ("fa ,  "Ip,  r)  - - — y  f  exp  [-  — — -4-1  I  |'/o|r/(l  +  r2)]  d/1 

°  2(1  +  r  )  1  2(1  +  r  )J 

- -  exp  I-  -ar-tl\  i  [V^r/(1  +  r2)]  d/3 

(1  +  r  )  1  2(1  +  r  )’  ° 


b/v  1+r 


(b/s/l  +  r2)  exp  j-  t-r.2>%+  ^H±  r¥j, 


IQ  [(ar /J  1  +  r2)  (b/V  1  +  r2)  J  d(b/x/l  +  r2) 


7  2  2 

1  a  r  +  b  ,  abr 

- — t-  exp - I  - y 

(1  +  r  )  2(1  +  r  )  °  1  +  r 
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Therefore, 


PQ(a,  b,  r) 


=  q( 


(.  ar  b  \ 

|a2r2  +  b2| 

1  !o| 

I  abr  1 

\l - 2  '  / - V 

1  2  exP_ 

(1  +  r  )  1 

‘  2(1  +  r2)  J 

‘1  +  r2 

(C-50) 


where  we  have  returned  to  the  parameters  a  =  sfo ,  b  =  •J'p  and  have  referred  to  Eq.  (C-4).  Max- 
imon's  result  [Eqs.  (C-2)  and  (C-3)]for  v  -  0  is  checked  by  substituting  the  second  series  of 
Eq.  (C-4)  in  Eq.(C-50). 
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Appendix  D 

COMMENTS  AND  ERRATA  FOR  APPENDIX  C 

Appendix  C  has  been  reproduced  in  its  entirety  from  Price  [  7  ] .  It  contains  the  neces¬ 
sary  derivation  to  complete  the  analysis  begun  in  the  present  report’s  text  of  probability  of 
error  with  no  smoothing  of  the  data. 

The  mathematics  in  Appendix  C  was  checked  and  found  correct  with  the  following 
exceptions: 

•  In  the  equation  which  follows  Eq.  (C-13)  the  exp  [-(a2  +  b2 )/2]  factor  in  the 
second  line  should  be  multiplying  the  entire  right-hand  side  of  the  equation  rather  than  just 
the  second  and  third  lines  as  shown. 

•  In  Eq.  (C-21)  the  exp  [-(a  +  0)]  factor  in  the  first  line  should  be  exp  [-(a  +  /3)/2] . 
In  neither  of  these  cases  does  the  error  propagate  to  the  equations  which  follow. 
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